arXiv:1508.01282vl [math.NA] 6 Aug 2015 


Approximating the Analytic Fourier Transform with the Discrete Fourier 

Transform 

Jeremy Axelrod* 

Department of Physics , University of California, Berkeley 
(Dated: 26 May 2015) 

The Fourier transform is approximated over a finite domain using a Riemann sum. This Rie- 
mann sum is then expressed in terms of the discrete Fourier transform, which allows the sum to be 
computed with a fast Fourier transform algorithm more rapidly than via a direct matrix multipli¬ 
cation. Advantages and limitations of using this method to approximate the Fourier transform are 
discussed, and prototypical MATLAB codes implementing the method are presented. 


I. Introduction 

The Fourier transform is a ubiquitous analytical 
mathematical tool. However, in many problems 
of interest analytic expressions for transformed 
functions do not exist, or the function to be trans¬ 
formed is only known at a set of discrete points 
as is the case for most real-world experimental 
data. In both of these cases, it is necessary to 
approximate the Fourier transform on a set of 
discrete points. This can be done by approximating 
the integral in the Fourier transform as a Riemann 
sum. Such a summation implemented as a single 
matrix multiplication by the vector of points to 
be transformed results in an undesirable algorithm 
complexity scaling of 0(N 2 ). Here, an algorithm 
is presented which allows this Riemann sum to be 
expressed in terms of the discrete Fourier transform 
(DFT), which can in turn be computed via a fast 
Fourier transform (FFT) with complexity scaling 
of O(NlogN) [1], The overall complexity of this 
algorithm also scales as 0(N log IV). 

II. Method 


A. Definitions 

Let the Fourier transform, /(w), of a function f{t') : 
R —» C be defined as 

(1) 

where 

(2) 

is the inverse Fourier transform of /, and a,b G 
R, b yf 0 are arbitrary constants chosen by conven¬ 
tion [2]. Let x := [xi, X 2 , ■ • ■ , £jv] be a vector of 
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length N, and let j, j', k, k' = 1, 2, 3,..., N. The dis¬ 
crete Fourier transform, X := [X\, X 2 , ••• , Xjv], 
of x is then defined to be such that 

N 

X k := £ a^e-W-1X*- 1 ) (3) 

r -1 


where 


N 
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fc' = 1 


is the inverse discrete Fourier transform of X. 


B. Approximating the forward transform 

Let the function f{t') be represented by a vector x 
of length N, so that Xj> = /(f',) where f'-, = r'(j' — 
1) +t[ constitutes a vector of evenly-spaced points 
beginning at t\. Then (1) can be approximated via 
a Riemann sum, F(ui), as 


/ M = F(u) 
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ibujt'.f 


( 4 ) 


Letting maxjf'-,}, min { tf ,} denote the maximum 
and minimum values in {if,}, respectively, then 

f(t')e ibut dt' 

when the limits of integration minji',} and 
maxji',} are held fixed so that \t’\ decreases as N 

increases. Therefore, F(ui) can be viewed as a Rie¬ 
mann approximation of the Fourier transform of the 
function 


lim F(lo) = 

TV—>oo 
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fo(f) := 



V t' € [min |i',} , max { t'y }] 
elsewhere 
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so that F(w) approximates /(w) well for large N, 
and for {i',} and f(t') such that 
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f(t )e zbut dt + / 


{*} 


f{t')e ibut 'dt' 
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ax{t',} 


f{t')e ibut 'dt' V i 


Computing the sum in (4) directly for TV values 
of w (denoted by w := [wi, u> 2 , •••■■, wjv]) can be 
written as a matrix multiplication: 


F(d3) 



gibtoit '2 
gibu)2t , 1 e ibuj2t' 2 

gibLOjyt^ gibco N^2 


e ibu>!t' N " 


Xi 

e ibbj 2 t' N 
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(5) 


Since multiplying a vector into an IV x N matrix in¬ 
volves N 2 multiplication operations and N(N — 1) 
addition operations, the asymptotic complexity of 
computing (5) scales as 0{N 2 ). However, “linearith- 
mic” 0(N log N) scaling can be achieved by express¬ 
ing F (w) in terms of the DFT as defined in (3) 
and then utilizing an FFT algorithm to compute the 
DFT [1], Referencing (3), in order to express F(u>) 
in terms of X, define 


so that 
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As such, F has been expressed in terms of the DFT, 
X, for oj = {uik} as defined by (6). However, it is 
sometimes desirable to retrieve the values of F for 
u £ {wfc}, e.g. u> < 0 if b < 0. To that end, note 
(6) that for m € Z, 
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so that 


F 


is periodic with period -A, with a 


phase shift of -F rnt\. Since oj spans an entire 
period (from 0 to F can be determined 

from F(u) for any integer multiple of 2 £ N . This 
is useful in practice, since it is often desirable to 


retrieve F on the interval u> £ [— \u nyq \, + |w ny9 |] 
for LJ nyq := — ^ the Nyquist frequency in order 
to display aliases in a more readily interpretable 
context-that is, negative frequency components will 
appear below u> = 0 instead of above u> = uj nyq . 
A comparison between the DFT and the Riemann 
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sum approximation to the Fourier transform is 
shown in figure 1. Clearly, if the vector being trans¬ 
formed is not entirely real, considering the DFT 
only below the Nyquist frequency index disregards 
possible asymmetries between positive and negative 
frequency components. Also, the vertical scaling of 
the Riemann sum approximation is independent of 
N unlike with the DFT, making direct comparisons 
of spectral power density between transforms of 
vectors of different lengths possible. 

C. Approximating the inverse transform 

The inverse transform can be treated in exactly the 
same way as the forward transform. Letting the 
function /(<*/) be represented by a vector X of length 
N so that X k > = f(eo' k ,) where co k , = W'(k' — 1) 
constitutes a vector of evenly-spaced points begin¬ 
ning at u)[, then (2) can be approximated via a Rie¬ 
mann sum, F(t), as 

m = m := \w\ x^e-^ (?) 

V ( 27r ) *;'=i , 


Then, after defining tj := -w^U - 1), 

F[tj) = \f^^ N \ W '\e~ ibU ' lti x j 

and for ( S Z, 

similarly to before. 

D. Invertibility of the approximate transforms 

It is of interest to know if the above Riemann sum 
approximations of the forward/inverse Fourier trans¬ 
forms are inverses of each other, i.e. is F applied to 
X = F(Hj), where F is applied to vector x. equal to 
x? Letting xy be defined on ty = t’( f—l)+t' 1: and 
using equations (4) and (7), the inverse transform of 
the forward transform of x can be written as 
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assuming that u)' k , = W'{k' — 1) + in accordance 
with equation (6) and the accompanying discussion 
on discrete shifts. The sum in parentheses in (8) can 
be evaluated as a geometric sum: 
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£(■ 

k'=l 


ibW'(t'., —t ) 


fc'-l 


1—e 


iNbW' (t'., -1 ) 


= ^ l—e 

N 


ibW'(t'.,-1 ) 


for e ibW ' (t i'- t) ± 1 
for j bw ' {t 'i'~ t] = 1 


Since l T 'l = N\W'\\b\ >. 

^iNbW' (ty — t) _ ^ 

^ibW'(tjf—t) _ ^sgn(bW , ')i2ivm/N ^ 

V m G {Z \ 0 | \m\ < N — 1} 

Therefore, 

V I r ibw\t'.,-t)\ k '^ 1 = fo for t’y ± t 
and so 


In order to consider invertibility, it must be that 
t € If this is the case, then V j' € 

{1,2,3,..., iV}, ty — t = m\r'\ for some m € 
{Z | \m\ < N — 1}. That is, the difference between 
t and ty is always an integer multiple of \t'\. Thus, 


F{t) = \t j \ \ W'\ Nxj> = xy when t = tf, 

which proves that the Riemann sum approximations 
are indeed inverses of each other. Since the method 
for calculating the Riemann sum approximation 
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FIG. 1. Comparison of the magnitude of the DFT (a) and Riemann sum approximation to the Fourier transform (b) 
for f(t) = sin(t) + 0.1e~ 2lt sampled at 201 evenly-spaced points on the interval [—100,100] with a = 0, b = —1. The 
DFT is more difficult to directly interpret because of aliasing across the Nyquist frequency at k = 101. 


of the inverse Fourier transform using the inverse 
DFT can only return t which are integer multiples 
of N w'b > transforms F(t) and F(ui) specified 
above are restricted to being each others’ inverse 
under the single condition that the input {t'.,} is of 
the form £'•, = r'(j' — 1 + Z) so that it is possible for 
the inverse transform to return t £ {/',}. 


III. Performance 


A prototypical MATLAB function, FT, which imple¬ 
ments the approximation method described above 
for the forward Fourier transform is given in Section 
Al. The Riemann sum approximation to the Fourier 
transform of /(£) = rect(i — 1) calculated using FT 
is shown in figure 2. The analytic Fourier transform 
of /(f) for a = 0, b = —1 is 


/M = 


1 


sm 


(l) 


V2tt 


;/2 


(9) 


and is also plotted in figure 2 for comparison. 
Referencing the code in FT, it is evident that adapt¬ 
ing the DFT to the Riemann sum approximation 
only incurs additional computational costs scaling 
as O(N) from the element-wise multiplications by 
phase factors. Since the FFT complexity scales as 
LT(TV log TV), the overall complexity of the adap¬ 
tation algorithm scales as O(NlogN). This is 
demonstrated in figure 3. 


IV. Conclusions 


It has been shown that the Riemann sum approxi¬ 
mation to the Fourier integral over a discrete finite 


domain can be expressed in terms of the discrete 
Fourier transform and is therefore calculable using a 
fast Fourier transform algorithm, which reduces the 
complexity of the problem from 0(N 2 ) for a direct 
matrix multiplication implementation of the sum to 
(H (TV log TV). The Riemann sum approximation is 
useful when a discrete approximation to the contin¬ 
uous Fourier transform is required, and it may be 
preferable to the discrete Fourier transform in some 
cases because it is directly interpretable in the con¬ 
text of the Fourier transform. The Riemann sum ap¬ 
proximation of the inverse Fourier transform applied 
to the Riemann sum approximation of the forward 
Fourier transform of a vector x returns the same vec¬ 
tor x —that is, the approximated transforms are still 
inverses of each other. However, the method pre- 



FIG. 2. The real (solid blue line) and imaginary (solid 
red line) parts of the analytic Fourier transform of f(t) 
given by equation (9) are plotted with the real (blue 
•) and imaginary (red •) parts of the Riemann sum ap¬ 
proximation as computed by MATLAB function FT. FT 
used 201 evenly-spaced samples of /(f) on the interval 
[—10,10] with a = 0, b = —1. The approximation per¬ 
forms well for |w| <C \ . 
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FIG. 3. Algorithm execution times as a function of input vector length for a) the Riemann approximation method 
as given in the function FT (see Section Al) and b) MATLAB’s built-in fast Fourier transform function, fft. c) The 
ratio of times in (a) to (b). The ratio approaches a constant value as N increases, indicating that the complexity 
of function FT scales similarly to fft -that is, as 0(N log N). The transforms were performed on the rectangular 
function f(t) = rect (t — 1), with sample points evenly-spaced on [—10,10]. 


sented here results in the domain of the transformed 
function being defined only on integer multiples of 
2^ where A is the spacing between values in the 
conjugate domain. Therefore, the transforms can 
only be implemented in such a way that they are 


each others’ inverse if the input domain is defined 
on integer multiples of some number. Future work 
could try to alleviate this integer multiple constraint, 
or could try express higher-order integral approxi¬ 
mations (like Simpson’s rule) of the Fourier integral 
in terms of the discrete Fourier transform. 
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Al. Appendix 

A simple implementation of the forward transform in MATLAB is given below. The function FT takes input 
vectors of length TV t, f and input scalars w_l, a, b (with w_l an integer multiple of t J^ n ) and returns output 
vectors of length N w,ff where w(l) = w_l. 

function [w,ff] = FT(t,f,w_l,a,b) 

N = numel(t); % size of input vectors 

tau = ((max(t)-min(t))./(N-l)); % input sampling period 

W = -2*pi/(tau*b*N) ; % define frequency spacing 

w_k = W.*(0:N-1); % define frequency vector beginning at zero 

X = sqrt(abs(b) ./((2*pi) ." (1-a)))*tau.*exp(li*b*w_k*t(1)) .*fft (f); 

% unshifted Riemann sum approximation using built-in "fft" function 

w = w_k+w_l; % shift frequency vector so that w(l) = w_l 

n = round(w_l/W); % number of indices that w is shifted from w_k by 

m = zeros(size(w_k)); % initialize vector of periodic shift indices 
if n > 0 

% determine m for each entry in w: 
idx = mod(n,N); 

m(N-idx:N) = floor(n./(N+l)) + 1; 
m(l:N-idx-l) = floor(n./(N+l)); 

inds = circshift(1:N,-idx,2); % circularly shift indices of X to 
% corresponding new w positions 
elseif n < 0 

% determine m for each entry in w: 
idx = mod(-n,N); 

m(l:idx) = -( floor(-n./(N+l)) + 1 ); 
m(idx+l:N) = - floor(-n./(N+l)); 

inds = circshift(1:N,idx,2); % circularly shift indices of X to 
% corresponding new w positions 

end 

phase = (exp (li* (2*pi/tau) *t (1) ) ) . ''m; % generate phases to shift by 

ff = phase.*X(inds); % phase shift elements of circularly-shifted X, output 

% answer 

end 
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The corresponding inverse transform function IFT is also given for convenience: 

function [t,f] = IFT (w, ff, t_l, a,b) 

N = numel(w); % size of input vectors 

W = ((max(w)-min(w))./(N-l)); % input sampling period 
tau = -2*pi/(W*b*N) ; % define time spacing 

t_j = tau.*(0:N-l); % define time vector beginning at zero 

x = sqrt(abs(b)./((2*pi).~(1+a)))*W.*exp(-li*b*t_j*w(l)).*N.*ifft(ff); 

% unshifted Riemann sum approximation using built-in "ifft" function 

t = t_j+t_l; % shift time vector so that t ( 1 ) = t_l 

n = round(t_l/tau); % number of indices that t is shifted from t_j by 

m = zeros(size(t_j)); % initialize vector of periodic shift indices 
if n > 0 

% determine m for each entry in w: 
idx = mod(n, N); 

m(N-idx:N) = floor(n./(N+l)) + 1; 
m(l:N-idx-l) = floor(n./(N+l)); 

inds = circshift(1:N,-idx,2); % circularly shift indices of X to 
% corresponding new w positions 
elseif n < 0 

% determine m for each entry in w: 
idx = mod(-n,N); 

m(l:idx) = -( floor(-n./(N+l)) + 1 ); 
m(idx+l:N) = - floor(-n./(N+l)); 

inds = circshift(1:N,idx,2); % circularly shift indices of X to 
% corresponding new w positions 

end 

phase = (exp(-li*(2*pi/W)*w(1))).~m; % generate phases to shift by 

f = phase.*x(inds); % phase shift elements of circularly-shifted x, output 

% answer 

end 


